Vertical Density Matrix Algorithm: A Higher-Dimensional Numerical 
Renormalization Scheme based on the Tensor Product State Ansatz 

Nobuya Maeshima, Yasuhiro Hieida^, Yasuhiro Akutsu, Tomotoshi Nishino^, and Kouichi Okunishi^ 

Department of Physics, Graduate School of Science, Osaka University, Toyonaka, 560-0043, Japan 
^Department of Physics, Graduate School of Science, Kobe University, Rokkoudai, 657-8501, Japan 
^Department of Physics, Faculty of Science, Niigata University, Igarashi, 950-2181, Japan 

(February 1, 2008) 

We present a new algorithm to calculate the thermodynamic quantities of three-dimensional (3D) 
classical statistical systems, based on the ideas of the tensor product state and the density matrix 
renormalization group. We represent the maximum-eigenvalue eigenstate of the transfer matrix 
as the product of local tensors which are iteratively optimized by the use of the "vertical density 
matrix" formed by cutting the system along the transfer direction. This algorithm, which we call 
vertical density matrix algorithm (VDMA), is successfully applied to the 3D Ising model. Using 
the Suzuki- Trotter transformation, we can also apply the VDMA to two-dimensional (2D) quantum 
systems, which we demonstrate for the 2D transverse field Ising model. 

PACS numbers; 05.10.-a ; 05.10.Cc ; 05.50.-Hq ; 64.60. Cn 



I. INTRODUCTION 

Since the density matrix renormalization group (DMRG) method was invented by S.R.White, [Q the method 
has been applied to various problems in one-dimensional (ID) quantum systems and two-dimensional (2D) classical 
systems. ||^ Such a great success of the DMRG has been stimulating us to extend the algorithm to the one which can 
handle higher-dimensional systems, principally 2D quantum systems and 3D classical systems. 

We should recall that, in the DMRG, the matrix-product structure of the wavefunction of the target states (usually 
the ground state or the maximum-eigenvalue eigenstate) is essential. Q From this point of view, the tensor product 
state (TPS) which is a natural higher-dimensional generalization of the matrix product state, should play a key 
role in the "higher-dimensional DMRG" . A simple but non-trivial example of the TPS is the ground state of the 
2D valence-bond-solid (VBS) type quantum spin systems, where the wavefunction is expressed as a product of local 
finite-dimensional tensors, with all the tensor indices being contracted. As for 3D classical statistical systems, 
the maximum-eigenvalue eigenstate of the layer-to-layer transfer matrix of the 3D classical system can exactly be 
represented as the TPS, if we allow the tensor dimension to be infinite. We should note that we can reduce the 
calculation of the expectation value of the TPS to a statistical average in a lower-dimensional classical system; 
{D -\- l)-dimcnsional classical (or D-dimensional quantum) problem reduces to a Z3-dimcnsional classical statistical 
problem. |^ In fact, the properties of the 2D VBS model have been studied in terms of a 2D vertex model associated 
with the TPS wavefunction. §,0 

When developing the TPS formulation, the most important step is the determination (optimization) of the local 
tensor in the TPS. In Refs. ||l^,|l^, direct variational formulations are employed for the optimization of the local 
tensors. As compared with these "direct" method, our novel algorithm given in the present paper for 3D classical 
system is more like the original DMRG. The local tensor is updated by the "block-spin basis transformation" along the 
vertical direction. Since this transformation is constructed in terms of the density matrix made along the "vertical 
direction" (transfer direction of the transfer matrix), we call this algorithm the vertical density matrix algorithm 
(VDMA). We apply the VDMA to the 3D Ising model and discuss its efficiency. We also report the application of the 
VDMA to the 2D transverse field Ising model, with help of the Suzuki- Trotter transformation p3[ |. 

This paper is organized as follows: In Sec.||, we briefiy explain the VDMA for 3D classical spin systems, taking the 
3D Ising model as an example. In Sec.pl| we show the numerical result for the 3D Ising model and the 2D transverse 
field Ising model. The last section is devoted to the conclusion. 

II. METHOD 

Let us consider the 3D Ising model on the simple cubic lattice of the size iV x iV x 2L in X, Y, and Z directions. 
Suppose that L and N are sufficiently large, and the neighboring Ising spins a and a' have ferromagnetic interaction 
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where i' — i + 1, j' = j + 1, and K — J/T. The locations of the spin variables are shown in Fig. |l|. 

For the book keeping, let us here introduce some notations for spin variables. Write the configuration of the four 
spins surrounding a plaquette in the XY plane as 



Wij} = [cTijCri'jai'j'atj'), 
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where the position of the plaquette can be labeled by the index ij. Then the Boltzmann weight is simply written as 
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Also for a spin layer in the XY plane, we denote the configuration of the N x N spins as, 
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Using these notations, the transfer matrix T from a layer [a] to the next layer [a] is written as: 



T[a\a]= n ^{J 



(4) 



In the product of Eq. (^, the spin variables are shared by the adjacent cubes in the diagonal direction and thus the 
Boltzmann weights make the checkerboard pattern in the XF-plane.(see Fig.l-(b)) 

Our goal is to evaluate the maximum-eigenvalue Amax of the transfer matrix T and the corresponding eigenvector 
iV'max), using the TPS representation of the eigenvector. In order to analyze the TPS structure of ji/'max), let us 
consider the power method briefly, which is the simplest but powerful technique to calculate Amax and IV'max)- Define 
the vector with 

IV'l) =T^|V'o), (5) 
. Then the maximum-eigenvalue eigenvector IV'max) is 



where |V'o) is an 'initial' vector that is not orthogonal to IV": 
obtained as 
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where the details of IV'o) are not important. 

In Fig. H, we show the graphical representation of |V'l), where we can see the structure of IVl) more clearly. As 
is seen in this figure, the L unit cubes are piling vertically up to the surface, and then the product of the vertically- 
arranged Boltzmann weights can be regarded as a local tensor; We define the local tensor at the ij-plaquette 
as, 
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where {cfj} is the spin configuration at the ij-plaquette on the "surface layer", and the auxiliary variable = 
i'^ij'^ij'^ij ■ ' ■ denotes the configuration for the spins under the surface. Using the local tensors defined above, 

IV'l) is represented as a TPS: 



(8) 
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Taking the limit L — > oo, we have the maximum-eigenvalue eigenvector IV'max), which is now represented as the 
product of the local tensor A^o ■ 
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From the practical view point of the usual power method, the eigenvector is improved with the relation ItpL^i) = 
T\ijjL) iteratively. In the framework of the TPS, this power method is reformulated as the recursion relation for the 
local tensor: 



with which we can carry out the iterative calculation, until gives a good approximation of Aoo- However it is 
generally difficult to store the tensor in the computer memory for a sufficiently large L, because if has M 
states, then the extended auxiliary variable ^^+^ has 2M states. 

In order to restrict the number of the auxiliary variable, we now import the idea of the DMRG into the TPS. The 
essence of the DMRG is that the increased number of states for ^^"""^ can be reduced, by using the "projection operator" 



generated from the "density matrix" . In the present case, the appropriate density matrix should be constructed for 
the spin variables i'^ij^^ i^ij'^^) i^i the vertical direction (we thus call this density matrix as "the vertical density 
matrix"). [Q Introducing the "transposed" local tensor: 
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the explicit form of the vertical density matrix is defined as (See Fig|^.), 
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where ^ denotes the configuration sum for the all spin variables except a^^i , ^^.^ , and Y[ means the product for the 
site index (i, j) except at (i,j) = (fc,Z), {k",l") with k" = k — 1 and I" = I — 1. In Eq. (pi]), we have also used the 
notation for the "checked spins" : 



(jk-i'' 

^k"l" 



Cfk"l" CTkl" (Jkl (^k"l 
ik"l" (,kl" Ckl ^k"l 



Further, we have omitted the sub(super)script L if it is apparent. 

For the case of the isotropic 3D Ising model, the vertical density matrix pki becomes independent of the site index 
kl in the thermodynamic limit. Thus we write pki simply as p hereafter. Moreover it should be noted that, for the 
isotropic case, the local tensors A and A satisfies the relation: 
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Following the spirit of the DMRG, we diagonalize p to have the eigenvalues Wc in the decreasing order wi > W2 > 
•(>0): 



(13) 



where U{a£,\^) is the eigenvector for w^. By taking U{a£,\^) with ^ e 1, • • • ,M, we construct the projection operator 
U, which is a 2M x M rectangular matrix. Operating U to the spin variables on each "edge" of A, we then make the 
renormalized local tensor A: 
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Using Eqs. (^ and (14) recursively, we can now calculate the effective local tensor Aoc- However we encounter 
another problem in this process; it is also a numerically heavy problem to compute the vertical density matrix with 
Eq. (P). 
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Let us next explain how to overcome the difhculty in calculating the vertical density matrix. The key idea is that 
we can consider Eq. (|l^) as a kind of 2D classical spin system with a point defect. To see it, we here define a new 
tensor 

'""^ '''' 

^ ^'^3 ) 

which is graphically represented in Fig. ^. 

Regarding the spin variable along the z axis as a 2M^ state single spin: 

(16) 

we can see that G{cry} becomes the Boltzmann weight for the 2D effective classical model. In fact, the vertical 
density matrix p can be expressed as the density matrix for the point defect in the 2D classical model: 




Pki 
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where the meaning of the prime at the summation and the product is the same as Eq. (|llD, and the "checked spin" 
is given by 

{d'k"i' 

We calculate each component of the vertical density matrix ([17|) as a statistical average of a special local observable 
in 2D. To this end we apply the corner transfer matrix renormalization group (CTMRG) JTsj l which has been known 
to be quite efficient for 2D classical statistical systems. 

Thus we have obtained a closed algorithm to calculate the local tensor A with Eqs. (^ and (p^, assisted by the 
CTMRG for the 2D effective classical model. We summarize the numerical procedure as follows: 

(a) For the free-boundary condition in the Z direction, define the initial local tensor Ai as 
For the Ferro boundary condition, define Ai as 

where {+1} means that the four Ising spins are aligned upward. 

(b) Define the effective Boltzmann weight of the 2D classical spin system with Eq. (p^). 

(c) Perform the CTMRG calculation for the 2D effective classical spin system with the Boltzmann weight G and 
obtain the vertical density matrix p. 

(d) Diagonalize p and construct the projection operator U. 

(e) Renormalize the local tensor A with Eq. (p^). 

(f) Return (b), until the local tensor A is converged. 

In this VDMA calculation, the accuracy is determined by the number of retained basis M for the auxiliary variables 
^ and 7], and m for the CTMRG calculation in 2D classical system. We can check the convergence of the computed 
quantities with respect to M and m. 
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III. RESULTS 



A. The 3D Ising model 

Fig. H shows the spontaneous magnetization (cr) calculated by using the VDMA. For comparison we also show the 
results of the 3D-version of the Kramers- Wannier (KW) approximation and the Talapov and Blote's Monte Carlo 
results For each Af, we have carried out the VDMA calculations for m = 4, 8, 12, and 16, where we can find 
the good convergence with respect to m. We observe that the convergence with respect to M is also sufficient in the 
off-critical region. Near the critical point, however, the magnetization becomes smaller as M is increased, implying 
that large M is needed for calculation in the critical region. 



B. The 2D transverse field Ising model 

Let us next consider the VDMA for the 2D transverse field Ising (TFI) model on the square-lattice, which is one 
of the fundamental quantum spin models in 2D. The Hamiltonian of the 2D TFI model is given by 

iJ = -J^afa|-r^af, (20) 

where J(> 0) is the ferromagnetic coupling constant, erf (a = x,y,OT z) are the Pauli matrices, and {i,j) denotes the 
nearest-neighbor pairs on the square lattice. The transverse field F induces the quantum fluctuation into the system. 
At zero temperature, the TFI model exhibits the quantum phase transition at F = Fc, below which the spontaneous 
magnetization (cr) behaves as (a) ^ (Fc — F)''. The critical field Fc and the critical exponent /3 have been estimated 
as Fc = 3.06 and (3 = 0.31 by the quantum Monte Carlo (QMC) simulation. |jl^ In the following we consider the 
VDMA for the 2D TFI model at zero temperature, and show the results of (cr). 

As was described in the previous section, the VDMA is formulated for the 3D classical systems. In order to apply 
the VDMA to the 2D TFI model, we map the model to the 3D anisotropic Ising model by using the Suzuki- Trotter 
transformation |l^. The partition function Z of the 2D TFI model is obtained as the limit of the 3D anisotropic 
Ising model: 



Z = lim Tr exp 

L — 'oo 



L L 

Kh ^ ^ cri,T<yj,T + Kv ^ ^ ai^rCi^T+i 

T=l (jj) T=l i 



(21) 



where ai^r is the Ising variable at the position i and imaginary time r. The effective couplings and in Eq. (|21 
are given by 



Kh = eJ, (22) 

2 ''^'^'^^ 



X„ = -ilog[tanh(eF)], (23) 



e = l/(Ti), (24) 

where the subscripts h and v denote the horizontal (X-Y) direction and the vertical (Trotter) one, respectively. We 
can perform the VDMA calculation for this anisotropic 3D Ising model. 

We should make a comment about the boundary condition, before proceeding to details. As was seen in the previous 
section, the open boundary condition is assumed in the VDMA. However, the periodic boundary condition is imposed 
along the Trotter direction in Eq. (^l|). As far as the zero-temperature properties are concerned, the boundary 
condition is inessential due to the double limit T ^ and L ^ 0, allowing us to apply the VDMA to Eq. (|2l]). |]l8|] 

For a fixed value of e, we calculate the magnetization (cr(e)) with the VDMA for the infinite volume. After obtaining 
(cr(e)) for various e, we take the e — > limit by extrapolation. In the actual calculation, we have observed the following 
e-dependence 

(cr(e)) = (cr(0)) + const x (25) 
which we adopt for the extrapolation formula. 
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In Fig. 1^, we show the (cr) .vs. plots at F = 3.0 as an example, which are obtained with the VDMA of the 
numbers of retained bases {M, m) = (2, 8) and (3, 18), where the convergence with respect to m is rapid. In the region 
F > 3.2, (cr) converges to smoothly. 

In Fig. 1^, we show the (cr)-F curve and the results of the series expansion for comparison In F < 2.6, the 

VDMA results are sufficiently reliable, where the good convergence about m and M can be seen. We can further 
see the good agreement with the series expansion in the small field region (F < 2.0). In the vicinity of the critical 
point, however, the calculated magnetization exhibits the M dependence. The roughly estimated critical field from 
the VDMA calculation is about 3.2, which is 4% larger than the QMC one. 



IV. CONCLUSION 



In this paper we have constructed a novel higher-dimensional numerical renormalization algorithm which utilizes 
the natural tensor-product form of maximum-eigenvalue eigenstate iV'max) of the transfer matrix. In our algorithm, 
called VDMA (vertical density matrix algorithm), the local tensor forming IV'max) is iteratively updated using the 
vertical density matrix. We have successfully applied the VDMA to 3D Ising model. The VDMA can also be applied 
to 2D quantum systems using the Suzuki- Trotter transformation to 3D classical statistical systems, which we have 
demonstrated for the 2D transverse field Ising model. Application of the VDMA to other 2D quantum systems such 
as the Heisenberg model is an important subject of study, which we are now undertaking. 
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FIG. 2. The graphical representation of the vector {tpL}- 
Below the y-plaquette, we find that the L unit cubes are 
piling up. 



FIG. 5. The spontaneous magnetization of the 3D Ising 
model. 
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FIG. 6. The extrapolation of the magnetization at F = 3.0. 
The solid lines are linear fits of data. 
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FIG. 7. e = Umit of {a)-r curve at T = 0. The arrow 
shows the critical field obtained by the QMC simulation [17]. 
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